Exterior of ellipsoid
Adapted from: (Floudas et al., 1999; Section 3.5) and (Lasserre, 2009; Table 5.1)
Introduction
Consider the polynomial optimization problem from (Floudas et al., 1999; Section 3.5)
A = [
0 0 1
0 -1 0
-2 1 -1
]
bz = [3, 0, -4] - [0, -1, -6]
y = [1.5, -0.5, -5]
using DynamicPolynomials
@polyvar x[1:3]
p = -2x[1] + x[2] - x[3]
using SumOfSquares
e = A * x - y
f = e'e - bz'bz / 4
K = @set sum(x) <= 4 && 3x[2] + x[3] <= 6 && f >= 0 && 0 <= x[1] && x[1] <= 2 && 0 <= x[2] && 0 <= x[3] && x[3] <= 3Basic semialgebraic Set defined by no equality
8 inequalities
4.0 - x[3] - x[2] - x[1] ≥ 0
6.0 - x[3] - 3.0*x[2] ≥ 0
24.0 - 13.0*x[3] + 9.0*x[2] - 20.0*x[1] + 2.0*x[3]^2 - 2.0*x[2]*x[3] + 2.0*x[2]^2 + 4.0*x[1]*x[3] - 4.0*x[1]*x[2] + 4.0*x[1]^2 ≥ 0
x[1] ≥ 0
2.0 - x[1] ≥ 0
x[2] ≥ 0
x[3] ≥ 0
3.0 - x[3] ≥ 0
We will now see how to find the optimal solution using Sum of Squares Programming. We first need to pick an SDP solver, see here for a list of the available choices.
import Clarabel
solver = Clarabel.OptimizerClarabel.MOIwrapper.OptimizerA Sum-of-Squares certificate that $p \ge \alpha$ over the domain S, ensures that $\alpha$ is a lower bound to the polynomial optimization problem. The following function searches for the largest lower bound and finds zero using the dth level of the hierarchy`.
function solve(d)
model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c, p >= α, domain = K, maxdegree = d)
optimize!(model)
println(solution_summary(model))
return model
endsolve (generic function with 1 method)The first level of the hierarchy gives a lower bound of -7`
model2 = solve(2)-------------------------------------------------------------
Clarabel.jl v0.6.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 9
constraints = 12
nnz(P) = 0
nnz(A) = 24
cones (total) = 2
: Zero = 1, numel = 4
: Nonnegative = 1, numel = 8
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 2.2646e+00 2.2646e+00 0.00e+00 5.04e-01 2.23e-01 1.00e+00 2.84e+00 ------
1 4.8146e+00 4.8773e+00 1.30e-02 1.11e-01 3.74e-02 2.33e-01 6.39e-01 8.00e-01
2 5.8894e+00 5.8986e+00 1.56e-03 1.43e-02 4.91e-03 3.33e-02 9.41e-02 8.57e-01
3 5.9970e+00 5.9972e+00 3.18e-05 2.67e-04 9.33e-05 6.56e-04 1.82e-03 9.81e-01
4 6.0000e+00 6.0000e+00 3.18e-07 2.67e-06 9.32e-07 6.56e-06 1.82e-05 9.90e-01
5 6.0000e+00 6.0000e+00 3.18e-09 2.67e-08 9.32e-09 6.56e-08 1.82e-07 9.90e-01
6 6.0000e+00 6.0000e+00 3.18e-11 2.67e-10 9.32e-11 6.56e-10 1.82e-09 9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 577μs
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -6.00000e+00
Dual objective value : -6.00000e+00
* Work counters
Solve time (sec) : 5.76540e-04
Barrier iterations : 6The second level improves the lower bound
model4 = solve(4)-------------------------------------------------------------
Clarabel.jl v0.6.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 82
constraints = 101
nnz(P) = 0
nnz(A) = 242
cones (total) = 10
: Zero = 1, numel = 20
: Nonnegative = 1, numel = 1
: PSDTriangle = 8, numel = (10,10,10,10,...,10)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 9.7795e-01 9.7795e-01 1.11e-16 6.53e-01 4.62e-01 1.00e+00 2.05e+00 ------
1 1.2254e+00 1.2350e+00 7.85e-03 1.33e-01 7.13e-02 1.75e-01 4.77e-01 7.96e-01
2 2.6674e+00 2.7439e+00 2.87e-02 6.88e-02 2.81e-02 1.71e-01 2.05e-01 7.36e-01
3 4.3661e+00 4.4172e+00 1.17e-02 1.60e-02 5.48e-03 7.83e-02 4.58e-02 8.30e-01
4 5.1102e+00 5.1418e+00 6.18e-03 7.99e-03 2.62e-03 4.73e-02 2.38e-02 5.82e-01
5 5.5386e+00 5.5472e+00 1.54e-03 2.33e-03 7.37e-04 1.33e-02 7.23e-03 7.49e-01
6 5.6457e+00 5.6476e+00 3.40e-04 5.12e-04 1.62e-04 3.00e-03 1.66e-03 8.07e-01
7 5.6906e+00 5.6907e+00 1.34e-05 1.98e-05 6.26e-06 1.19e-04 6.47e-05 9.72e-01
8 5.6922e+00 5.6922e+00 8.64e-07 1.27e-06 4.02e-07 7.63e-06 4.16e-06 9.40e-01
9 5.6923e+00 5.6923e+00 1.21e-08 1.76e-08 5.58e-09 1.06e-07 5.77e-08 9.87e-01
10 5.6923e+00 5.6923e+00 2.66e-10 3.81e-10 1.21e-10 2.33e-09 1.25e-09 9.78e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 4.18ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -5.69231e+00
Dual objective value : -5.69231e+00
* Work counters
Solve time (sec) : 4.17628e-03
Barrier iterations : 10The third level improves it even further
model6 = solve(6)-------------------------------------------------------------
Clarabel.jl v0.6.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 451
constraints = 506
nnz(P) = 0
nnz(A) = 1376
cones (total) = 10
: Zero = 1, numel = 56
: PSDTriangle = 9, numel = (55,55,10,55,...,55)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 6.2703e-01 6.2703e-01 0.00e+00 7.56e-01 5.36e-01 1.00e+00 1.48e+00 ------
1 7.4953e-01 7.6135e-01 1.18e-02 1.58e-01 9.97e-02 2.39e-01 4.26e-01 8.06e-01
2 1.2333e+00 1.3121e+00 6.39e-02 4.31e-02 2.64e-02 1.73e-01 1.70e-01 8.68e-01
3 1.9671e+00 2.0062e+00 1.99e-02 1.22e-02 4.68e-03 6.32e-02 3.87e-02 8.18e-01
4 2.7047e+00 2.7536e+00 1.81e-02 8.45e-03 2.26e-03 6.99e-02 1.94e-02 7.43e-01
5 3.0495e+00 3.0614e+00 3.90e-03 3.75e-03 9.47e-04 2.13e-02 8.43e-03 9.22e-01
6 3.4706e+00 3.4774e+00 1.96e-03 9.59e-04 2.24e-04 9.63e-03 1.95e-03 8.04e-01
7 3.7081e+00 3.7129e+00 1.30e-03 4.91e-04 1.15e-04 6.51e-03 9.26e-04 6.34e-01
8 3.9158e+00 3.9173e+00 3.71e-04 1.29e-04 3.04e-05 1.93e-03 2.38e-04 8.24e-01
9 3.9469e+00 3.9488e+00 4.92e-04 8.66e-05 1.55e-05 2.34e-03 1.30e-04 6.86e-01
10 4.0073e+00 4.0086e+00 3.31e-04 3.98e-05 5.09e-06 1.52e-03 4.86e-05 7.76e-01
11 4.0267e+00 4.0275e+00 1.76e-04 1.93e-05 2.66e-06 8.15e-04 2.55e-05 7.35e-01
12 4.0582e+00 4.0584e+00 3.91e-05 3.52e-06 4.53e-07 1.79e-04 4.47e-06 9.70e-01
13 4.0660e+00 4.0660e+00 1.14e-05 9.59e-07 1.23e-07 5.20e-05 1.22e-06 8.03e-01
14 4.0671e+00 4.0671e+00 6.05e-06 4.79e-07 6.23e-08 2.75e-05 6.16e-07 6.69e-01
15 4.0683e+00 4.0683e+00 7.91e-07 6.17e-08 8.06e-09 3.60e-06 7.98e-08 8.78e-01
16 4.0684e+00 4.0684e+00 3.27e-07 2.44e-08 3.20e-09 1.48e-06 3.17e-08 7.67e-01
17 4.0685e+00 4.0685e+00 6.90e-09 5.14e-10 6.74e-11 3.12e-08 6.67e-10 9.80e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 47.1ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -4.06848e+00
Dual objective value : -4.06848e+00
* Work counters
Solve time (sec) : 4.71290e-02
Barrier iterations : 17The fourth level finds the optimal objective value as lower bound.
model8 = solve(8)-------------------------------------------------------------
Clarabel.jl v0.6.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 1813
constraints = 1938
nnz(P) = 0
nnz(A) = 5689
cones (total) = 10
: Zero = 1, numel = 126
: PSDTriangle = 9, numel = (210,210,66,210,...,276)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 5.6617e-01 5.6617e-01 4.44e-16 8.17e-01 7.34e-01 1.00e+00 1.32e+00 ------
1 6.4674e-01 6.4734e-01 6.01e-04 1.62e-01 1.39e-01 2.54e-01 4.12e-01 7.90e-01
2 8.6375e-01 9.4985e-01 8.61e-02 6.23e-02 5.34e-02 2.30e-01 2.32e-01 7.49e-01
3 1.3638e+00 1.3968e+00 2.41e-02 1.40e-02 8.73e-03 6.36e-02 5.21e-02 9.26e-01
4 2.0678e+00 2.0994e+00 1.53e-02 5.77e-03 2.46e-03 4.79e-02 1.82e-02 7.97e-01
5 2.3355e+00 2.3550e+00 8.36e-03 3.83e-03 1.37e-03 3.02e-02 1.09e-02 5.39e-01
6 2.7753e+00 2.7885e+00 4.75e-03 1.63e-03 4.62e-04 1.85e-02 4.07e-03 7.19e-01
7 3.1431e+00 3.1447e+00 5.16e-04 5.25e-04 1.36e-04 3.42e-03 1.30e-03 9.52e-01
8 3.2519e+00 3.2544e+00 7.58e-04 3.30e-04 8.17e-05 3.94e-03 7.58e-04 6.30e-01
9 3.3335e+00 3.3356e+00 6.53e-04 2.40e-04 5.77e-05 3.31e-03 5.40e-04 5.27e-01
10 3.5632e+00 3.5650e+00 5.04e-04 7.97e-05 1.71e-05 2.29e-03 1.62e-04 7.81e-01
11 3.6267e+00 3.6288e+00 5.82e-04 6.38e-05 1.32e-05 2.58e-03 1.22e-04 4.53e-01
12 3.8605e+00 3.8614e+00 2.08e-04 2.21e-05 4.24e-06 9.82e-04 4.09e-05 7.15e-01
13 3.9773e+00 3.9775e+00 2.96e-05 3.42e-06 6.21e-07 1.45e-04 6.28e-06 8.75e-01
14 3.9966e+00 3.9967e+00 5.04e-06 4.70e-07 8.32e-08 2.39e-05 8.76e-07 9.42e-01
15 3.9996e+00 3.9996e+00 6.22e-07 5.44e-08 9.53e-09 2.92e-06 1.02e-07 9.37e-01
16 4.0000e+00 4.0000e+00 5.82e-08 5.04e-09 8.83e-10 2.73e-07 9.46e-09 9.08e-01
17 4.0000e+00 4.0000e+00 1.63e-08 4.20e-09 2.64e-10 7.72e-08 2.73e-09 8.28e-01
18 4.0000e+00 4.0000e+00 2.64e-09 6.75e-10 4.27e-11 1.25e-08 4.42e-10 8.43e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 901ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -4.00000e+00
Dual objective value : -4.00000e+00
* Work counters
Solve time (sec) : 9.01331e-01
Barrier iterations : 18This page was generated using Literate.jl.